NCCS -5^/ 1 


M 


L ' " >V 


Summary of Work for 
Joint Research Interchanges 
with DARWIN Integrated Product Team 

1998 


Lambertus Hesselink 
Stanford University 
Stanford, California 94305-4035 
bert@kaos.stanford.edu 
Phone: 650-723-4850 Fax: 650-725-3377 

March 1999 


1 Introduction 

The intent of Stanford University’s SciVis group is to develop technologies that 
enabled comparative analysis and visualization techniques for simulated and 
experimental flow fields. These techniques would then be made available un- 
der the Joint Reasearch Interchange for potential injection into the DARWIN 1 
Workspace Environment (DWE). In the past, we have focused on techniques 
that exploited feature based comparisons such as shock and vortex extractions. 
Our current research effort focuses on finding a quantitative comparison of gen- 
eral vector fields based on topological features. Since the method relies on 
topological information, grid matching and vector alignment is not needed in 
the comparison. This is often a problem with many data comparison techniques. 
In addition, since only topology based information is stored and compared for 
each field, there is a significant compression of information that enables large 
databases to be quickly searched. This report will briefly (1) describe current 
technologies in the area of comparison techniques, (2) will describe the theory 
of our new method and finally (3) summarize a few of the results. 


2 Comparison Techniques 

There exist a variety of comparison techniques for vector fields. These tech- 
niques basically fall into three general categories: Image, data, and feature 
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extraction based comparisons. In most of these cases, comparisons are made 
visually [1]. Image based comparisons work on the computer generated image. 
Often times, a numerical data set is converted into an image that simulates an 
experimental visualization technique (computational flow imaging). This may 
be easier than extracting a vector field from an image, such as Schlieren. How- 
ever, visualizing a field in 3-D is quite difficult. Often times, these techniques 
are limited to two dimensions. In addition to side-by-side comparison of images, 
other techniques include image fusion, and Fourier analysis [2]. 

Data level comparison techniques operate directly on the raw data. An accu- 
rate comparison requires proper grid alignment which can involve problematic 
interpolation between two fields [3]. 

The last comparison category is the extraction of features. Typically fea- 
tures are flow specific such as vortex cores, shock surfaces, or topology. Often 
times there is a geometric representation of the feature and possibly a semantic 
representation of the system which can be compared using a pattern recogni- 
tion technique [4]. This may lead to more robust comparisons. Past study in 
our group has focused on the geometric structure of vector fields [5]. However, 
this geometric structure can be visually deceiving since two vector fields may 
have the same underlying topological structure but are dissimilar in appear- 
ance [6]. Therefore, a quantitative measurement for comparison of vector fields 
is essential. 


3 Description of a Vector Field 

A 2-D vector field can be described as a system of two simultaneous differential 
equations having the following form: 

v x = ^ = F (x,y) ( 1 ) 

Vy = ^=G{x,y) 

where F and G are continuous and have continuous partial derivatives in some 
region D. 

A vector field is typically described by the number, type, and arrangement 
of critical points (or equilibrium points). These points are where the system is 
defined to be F(x, y) = 0, G(x, y) = 0. The number and nature of critical points 
will not change under continuous transformation. A critical point is said to be 
isolated or simple if there is an open neighborhood around it that contains no 
other critical points. For this report, we focus entirely on simple critical points. 
The global topology of the vector field is defined as the critical points and the 
set of their connecting streamlines. These streamlines (separatrices) divide the 
field into regions that are topologically equivalent to uniform flow. Hence, only 
the topology is needed to reconstruct the field and therefore is useful as a means 
of differentiating vector fields. 
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Figure 1: basic patterns for simple critical points 


3,1 Classification of Simple Critical Points 

The behavior of the flow about a critical point can be analyzed by investigating 
the streamlines in the neighborhood of the critical point. If we are sufficiently 
close to the critical point (say a distance dx,dy away) in most cases a first order 
Taylor series expansion of the velocity field is sufficient: 

v x (dx,dy) ss ^ dx + ~^ d V ( 2 ) 

v v (dx,dy) « ~^ dx + 

Hence, the flow pattern is completely determined by the 2X2 Jacobian ma- 
trix, Jij = 1^*- (i, j = 1,2) evaluated at the critical point location. The various 
patterns formed in the phase-plane space can be seen by analyzing the eigenval- 
ues of the Jacobian. The patterns are sketched in Figure 1. Notice a positive or 
negative real part (denoted by a) is indicative of repelling/attracting behavior. 
And if an eigenvalue has an imaginary part (/? < 0), it indicates circulation 
about the point, otherwise asymptotic behavior is exhibited. 


4 Vector Field Representation 
using Clifford Algebra 

In [7] [8], Sheuermann et al. introduced Clifford algebra for vector field visu- 
alization. Clifford algebra provides a nice way to describe the relation between 
real and complex numbers in 2D space. The vector fields are defined over a 
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complex field in this algebra and the nonlinear vector fields are represented as 
multiplications of linear fields. 

For the Euclidean plane we get a 4-dimensional R-algebra G 2 with the basis 
l,ei,e 2 ,i = eie -2 as a real vector space. Multiplication is defined as associative, 
bilinear and by the equations 
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4.1 Vector fields in Clifford space 

A Clifford vector field is just a multi vector field with values in R 2 C G 2 

v : R 2 -> R 2 C G 2 (H) 

Let z = x + iy, z = x — iy be complex numbers in the Clifford algebra. This 
means 

x=-(z + z) (l 2 ) 

y = ^( z ~z) ( 13 ) 

We get 

v(r) = vi{x,y)ei +v 2 (x,y)e- 2 

= Vl G (2 + " ) 4 ( "“" ) ) ei 

~iv 2 Q (z + *) > Yi ^ _ ei 
= E(z,z)e 1 (14) 
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Table 1: Classification of Critical Points using a 0 values 


a 

0 

Type 

a 

0 

Type 

= 0 

< 0 

Center 

— 

>0 

Saddle \0\ > |a| 

> 0 

< 0 

Repelling Focus 

< 0 

< 0 

Attracting Focus 

> 0 

= 0 

Repelling Star 

< 0 

= 0 

Attracting Star 

> 0 

> 0 

Repelling Node \0\ < |a| 

< 0 

>0 

Attracting Node \0\ < \a\ 


Generally, a linear vector field can easily be shown as: 
v(r) = E(z,z)e i 

= (az + bz + c)ei (15) 


where a, 6 ,c C (7. 

Let £ : C 2 — > C be the polynomial so that v = £(z,z)ei. Let F* : C 2 —> 
C, k = 1, . . . , n be the irreducible components of E , so that F(z, z) = n*=i 
then an arbitrary polynomial vector field with isolated critical points can be 
expressed as: 

v(r) = E(z y z)ei 

n 

= nc.« + b k z + c k )e i (16) 

k = l 

where z* is the unique zero of a*z + 6 *z -f c*. 

5 a-/3 Space and its Use as a Metric 

For a linear vector field v — (az + bz T c)ei, let a = a\ + a^i and b — b\ + b^i- 
Eigenvalues of the Jacobian around its critical point zo are Ai = b\ T \J |a |‘ 2 — 6 2 
and A ‘2 = b\ — \/|a | 2 — 6 2 - 

Let a = b\ and 0 — sign(\a\ 2 — 6 2 )\/[H^ — 6 * 2 1, criteria for basic patterns of 
simple critical points are: 

Selection of a and 0 as shown in Figure 1 and delineated below can be 
mapped to 01 , 61 , 02,62 to yield any desired field: 

Notice our definition of saddle is more relaxed than shown in Figure 1. The 
values of a and 0 determine the type of critical point but it is not sufficient to 
be used as a metric to differentiate between two types of critical points. So we 
introduce a new a-0 space where the 8 simple critical points are mapped onto 
the a,0 axes at their respective (a,/?) points. Vectors in this space obey all 
the rules defined for a regular 2-D Euclidean space. All points in this space are 
normalized as follows: 


a = 


s/a?+g? 


& 


-J— 

y/a*+0 2 


(17) 


It is shown in [9] that the actual values of a and 0 do not determine the por- 
trait of the critical point only the ratio between them. Hence, this normalization 
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Figure 2: Basic patterns for critical points in a-0 space; C for center, RN for 
node, AN for attracting node, RF for repelling focus, AF for attracting focus, 
St for star, Sa for Saddle and R for regular point 


maps all points onto a unit circle (Figure 2) and thereby provides a means of 
relatively quantifying the difference between various points. Also note that a 
regular vector field with no critical points, v = const • ei has a = 0 and 0 = 0 
and sits at the origin of the unit circle. For the remainder of the report, a and 
(3 values will be assumed normalized. 

A multiple point with a set of a’s and 0 } s corresponds to a set of points in 
the q — 0 space. For example, v = z 2 e i is a dipole which has two (1,0) point 
in a - 0 space; and v = (z - (2 + 2 i)z 4- c\){z + (2 + 2 i)z + c 2 )e i has one point 
at (— and another point at in a — 0 space. 

6 Earth Mover’s Distance 

6-1 EMD analysis 

The Earth Mover’s Distance is first introduced in [10] [11] for content-based 
image retrieval in a large data base. It is used to compute the minimal amount 
of work that must be performed to transform one feature distribution into the 
other. Feature distribution in [10] [11] are the color and texture signatures of 
an image. 

After careful study, we found that the EMD concept can be used to compute 
the differences between vector fields. Here, the feature distribution is redefined 
as the characteristics of a vector field. 

Definition 1 (feature distribution) A feature distribution for a vector field 
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is the set of a and (3 values associated with the vector field's critical points: 

{(^ 1 j (3\ (^2 •> ^ 2 ) j • • • 5 (^n ? fin) } • 

Definition 2 (Energy) The energy for a vector field is: 

n 

Energy = , ^(ajf +/??), 

\ i= 1 

where n is the total number of critical points in this field. 

This energy here is a quantity that characterizes the critical points of a vector 
field. It is different from the physical energy. The concept "work” is used to 
measure the energy differences between two vector fields or the amount of energy 
used to transform one vector field into the other. 

Definition 3 (Work) For two vector fields with feature distributions 

{ (ck 1 , (3 \ ) , (02 j (3 2 )?-••) (<*n ? (3n ) } 

and 

{(a;,^),(af 2 ,^),...,«,/3;)}. 

The amount of work nec essary for transforming one ve ctor field into the other 
is defined as: Work = ” a I) 2 + (Pi “ $) 2 )* 

Intuitively, given two feature distributions, one distribution can be seen as a 
set of discrete point-objects with a certain amount of mass of earth spread in 
space, the other as a collection of holes in the same space. The work measures 
the least amount of energy needed to fill the holes with earth and is called the 
Earth Mover’s Distance (EMD). Computing the EMD is based on a solution to 
the old transportation problem from linear optimization [12]. This is a bipartite 
network flow problem which can be formalized as the following linear program- 
ming problem: Let I be a set of suppliers, J a set of consumers, C{j the cost to 
ship a unit of supply from i € / to j E J 



and it is the same as the Euclidean distance dij = ||vj — t?j|| in a-(3 space. A 
critical point either exists as a whole or does not exist, it can not be split. In 
this case, the transportation problem has the property that the optimal flow fij 
can only be 0 or 1 [13]. We want to seek a set of fij that minimizes the overall 
cost: 

EM D(x,y) — m.in^2 Cijfij (18) 

iel jtJ 
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subject to the following constraints: 


fij > o ie i, j e J 

(19) 

/ , fij = Vii j ^ J 

(20) 

i€l 


^ ] fij = x ii i £ I 

(21) 

j€J 


^2yj = 

(22) 

jeJ iel 



Where Xi is the total supply of supplier i and yj is the total capacity of consumer 
j. Constraint (19) allows shipping of supplies from a supplier to a consumer 
and not vice versa. Constraint (20) forces the consumers to fill up all of their 
capacities and constraint (21) limits the supply that a supplier can send as a 
total amount. Constraint ( 22) is a feasibility condition that ensures that the 
total demand equals the total supply, in other words, the distributions have the 
same overall mass and the EMD is a true metric [10]. 

It is likely that a set of vector fields will not have the same number of 
distributions. In order to satisfy constraint (22), we can create regular points 
to make the supply equal the demand without changing the vector fields. For 
example, if the supplier field contains 3 critical points 

3 

v = JJ(a*£ -F M 4- a)e i (23) 

i = 1 

and the consumer field contains 5 critical points 

5 

l>' = Y[( a 'i z + b 'iZ + c 'i) e l ( 24 ) 

j = i 

The supplier side has two fewer points in the a — (3 space. Now let 

3 

v = J"J(a;z -F ha -F c*) • 1 • lei (25) 

i = 1 

and the vector field remains unchanged. However, now we have two more regular 
points corresponding to 1 with a = 0 and (3 = 0, and both the supplier and 
the consumer have 5 points in their feature distributions. All the conditions are 
satisfied, and we are ready to compute the EMD for these two fields and find 
out the dissimilarity between them. 

In order to evaluate the meaningfulness of our new metric, we use Multidi- 
mensional MDS) [14] [15] to embed the vector fields in a two-dimensional 

Euclidean space so that distances in the embedding are as close as possible to the 
true EMDs between vector fields. The MDS is introduced in the next section. 
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Figure 3: EMD of the flow over an a) airfoil b)circular cylinder. 


7 Display of EMDs for a Large Set of 
Vector Fields 

The above discussions are for comparison of a pair of vector fields. If there 
exist a large set of vector fields and we want to compare their topologies, it is 
necessary to display them in a more meaningful way than a ID list sorted by 
their EMDs. Yossi Rubner et al. have used Multidimensional Scaling Method 
(MDS) [10, 11] to display a set of images on a 2D map. Given n objects in a 
high dimension, the MDS method computes a configuration in a lower dimension 
space such that the distance between every pair of objects in this low dimension 
space best matches the real distance in the high dimension. Inspired by their 
work, we compute the EMDs between every pair of vector fields and position 
the vector fields on a 2D map such that the distances between the vector fields 
match their EMD values as close as possible. 

8 Application: Flow over an Airfoil and Cylin- 
der 

Rogers and Kwak computed the flow past a 2-D airfoil at —90° angle of at- 
tack [16]. The model was of interest since the flow of the wake of an XV- 15 Tilt 
Rotor aircraft degraded the lifting capability during hover. An incompressible, 
time accurate, Navier-Stokes code with artificial compressibility at a Reynolds 
number of 200 was used to compute the flow over a NACA 64A223M airfoil. 
Fifty frames were computed. During this time the flow entered into a peri- 
odic vortex shedding cycle. Earth mover’s distance was computed over the 50 
frames. The plot in Figure 3a depicts the EMD comparison of frame 1 with 
the remaining 49 frames. At frame 1, the EMD is zero, which is expected since 
the work required to convert a frame into itself is zero. The periodic nature 
is apparent. We see a repetition approximately every 17 frames. Also we see 
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a sudden EMD rise when comparing frame 1 with frame 8 indicating a signif- 
icant topological feature difference. Frame 1 contains three critical points: an 
attracting/repelling focus and a saddle. Frame 8 contains 5 critical points: two 
saddles, an attracting/repelling focus pair and a node. Since the flow is incom- 
pressible, the velocity divergence, V • ff, is expected to be zero everywhere in 
the flow. Hence only saddles and centers are to be extracted. However, it is 
common due to numerical computation for a to not exactly be zero, however, 
we should not expect to find a node. Upon closer examination of the data, the 
velocity divergence for certain frames is not zero near the tips of the foil where 
nodes are being extracted. We believe that the flow solver may not have fully 
converged and therefore we see this sudden jump of discontinuity in the field. 
The LIC images of frames 1, 34, and 50 are depicted in Figure 4 have very 
similar earth mover’s distance and as can be seen look nearly identical. Frame 

8 differs from the others due to its variation in topology (formation of nodes) 
and is apparent in the figure. 

We contrast the flow over an airfoil with the flow over a circular cylinder 
simulated by Rogers and Kwak under the same flow conditions [16]. In this case, 
the flow is divergence free and the EMD values are quite similar. Thirty frames 
were computed capturing a complete cycle of vortex shedding. As can be seen 
from the plot in figure 3b, Frames 1 and 16 have nearly identical EMD values 
leading one to believe the period to be every 15 frames. Due to the symmetry of 
the flow, this is not far from the truth. In fact the flow produces a mirror image 
of itself every 15 frames as it sheds the alternate vortex and hence leads to the 
same topology. Figure 5b depicts the alternate vortex being shed to the image 
found in figure 5a. Furthermore, we see from figure 3b an increase in EMD 
value for frame 3. This increase is due to the dissipation of the saddle-center 
pair as it moves dow r n stream (figure 5c). The EMD drops by frame 6 as the 
next saddle-center pair is shed (figure 5d). By frame 16, the saddle-center has 
moved down stream such that the a, /? values are nearly identical to frame 1. 
Two frames later the saddle-center dissipate and the cycle repeats. 

9 Discussion 

We have demonstrated the effectiveness of topology based feature comparisons 
for vector fields. The use of a quantitative measure between fields provides the 
means for fast automated comparisons as well as an indepth study of flow fields 
as demonstrated with the time history data. For the airfoil data, we have shown 
the effectiveness of the method as a diagnostic tool. The clear EMD difference 
provides an immediate alert into calculation problems for particular frames. For 
the cylindrical data, the periodic nature of the flow was revealed. The EMD 
difference also provides insight into the evolution of the flow field. 

We currently are researching techniques to extend this method to three di- 
mensional vector fields and eventually to tensor fields. 
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Figure 4: Topologically similar frames 1, 34 and 50 of flow about an airfoil. 
Frame 8 is topologically dissimilar. 
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Figure 5: Topologically similar frames 1 and 16. Frame 4 depicts the down 
stream dissipation of saddle-center pair producing a larger EMD. Frame 6 de- 
picts formation of new saddle-center pair at the bottom of the cylinder. 
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